MVP 1

Advanced Geocomputation

Author

Matt Braaksma

Description

imports

Design Framework

  • Problem: Write a 2-3 sentence description of the problem you are solving.
    • OpenStreetMap provides layer of land use/land cover data. How does this compare to land use/land cover data derived from remote sensing data?
  • Solution: Write a 2-3 sentence description of the robust solution that would solve the problem. Note, this is not the MVP solution. This is the complete solution if you had unlimited resources.
    • Answering this question requires accessing both OSM LULC data and a remote sensing data set. Next a correspondence is required to ensure the proper classes are being compared. Then spatial correlation analysis can be conducted. Finally, other data could be included to determine whether the correlations are driving by other socioeconomic or geographic determinants.
  • Challenge: Write a 1-2 sentence description of the top-most challenge to create your solution.
    • Sourcing the data and ensuring that the different pieces work well together.
  • Spec list: Brainstorm all of the features/capabilities your solution will need. Rate each feature/capability on their value (how useful is it?) and effort (how much time to create it?). List your top 7 features and their rated value and effort as High, Medium, Low (H, M, L) for your spec.
Spec Value (H, M, L) Effort (H, M, L)
Get OSM data H L
Get remote sensing data H M
Create class correspondence M M
Rasterize OSM data H H
Compare OSM and remote sensing rasters H M

Review the spec list. Cross of features that are too high effort for too low value until you have a list of features that can be implemented in 2 weeks. That is your MVP! Success: Now that you have your MVP spec. What is the key metric for success? If you achieve XYZ you know that your MVP was successful. Write 1 sentence outlining your 1 key metric. Reflection: After you implement your MVP, write 6-7 sentences reflecting on your experiences in developing and running your MVP. Was it successful? What did you learn? Would you have done something differently if you were to repeat the project?

Code

Import modules

Code 
from pyrosm import OSM
from pyrosm import get_data
import matplotlib.pyplot as plt

import osmnx as ox
import geopandas as gpd
from osgeo import gdal
import pygeoprocessing as pygeo
Code 
# Get the geographic boundary for St. Louis County, Minnesota
gdf = ox.geocode_to_gdf("Lakeville, Minnesota, USA")
gdf_path = '/Users/mbraaksma/Files/advgeocomp2024/advgeocomp2024/mvp01/data/lakeville.gpkg'
# albers_crs = {
#     'proj': 'aea', 'lat_1': 29.5, 'lat_2': 45.5, 'lat_0': 23, 
#     'lon_0': -96, 'x_0': 0, 'y_0': 0, 'datum': 'WGS84', 'units': 'm', 'no_defs': True
# }
# gdf = gdf.to_crs(albers_crs)
gdf.to_file(gdf_path, driver='GPKG')
Code 
# File paths
input_raster = '/Volumes/T7 Touch/Research/Reservation_Land/NLCD/raw/nlcd_2021_land_cover_l48_20230630/nlcd_2021_land_cover_l48_20230630.img'
input_vector = '/Users/mbraaksma/Files/advgeocomp2024/advgeocomp2024/mvp01/data/lakeville.gpkg'
output_raster = '/Users/mbraaksma/Files/base_data/lulc/nlcd/nlcd_clipped.tif'

# Open the raster
raster = gdal.Open(input_raster)

# Use gdal.Warp to clip raster with the vector file
gdal.Warp(
    output_raster,
    raster,
    format='GTiff',
    cutlineDSName=input_vector,  # Path to vector file
    cropToCutline=True,          # Only crop the raster to the extent of the vector
    dstNodata=-9999,             # Set nodata value for the output raster
)

# Close the dataset
raster = None
Warning 1: for band 1, destination nodata value has been clamped to 0, the original value being out of range.

Import OSM LULC #2

Code 
gdf = ox.features_from_place("Lakeville, Minnesota, USA", tags = {"landuse": True})
gdf.explore(column='landuse')
/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/predicates.py:798: RuntimeWarning:

invalid value encountered in intersects

/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/set_operations.py:340: RuntimeWarning:

invalid value encountered in union

/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/predicates.py:798: RuntimeWarning:

invalid value encountered in intersects

/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/predicates.py:798: RuntimeWarning:

invalid value encountered in intersects

/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/predicates.py:798: RuntimeWarning:

invalid value encountered in intersects

/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/predicates.py:798: RuntimeWarning:

invalid value encountered in intersects
Make this Notebook Trusted to load map: File -> Trust Notebook

Import ESA LULC

Code 
lulc_path = '/Users/mbraaksma/Files/base_data/lulc/nlcd/nlcd_clipped.tif'
lulc_data = gdal.Open(lulc_path)

lulc_array = lulc_data.GetRasterBand(1).ReadAsArray()

# Close the dataset
lulc_data = None

# Plot the raster data
plt.figure(figsize=(10, 8))
plt.imshow(lulc_array, cmap='viridis')  # Adjust colormap if needed
plt.colorbar(label='Raster Values')
plt.title('Raster Plot')
plt.show()